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Abstract 

In this paper, the effect that produces the local linearization of the embedded 
Runge- Kutta formulas of Dormand and Prince is studied. For this, embedded 
Locally Linearized Runge-Kutta formulas are defined and their performance is 
analyzed by means of numerical simulations. The results show that the locally lin- 
earized formulas exhibit significant higher accuracy than the original ones, which 
implies a substantial reduction of the number of time steps and, consequently, a 
sensitive reduction of the overall computation cost of their adaptive implementa- 
tion. 

Keywords: Differential equation; Local Linearization; Runge-Kutta; Numerical 
integrator 



1. Introduction 

It is well known (see, e.g., [3, 21, 19]) that conventional numerical schemes 
such as Runge-Kutta, Adams-Bashforth, predictor-corrector and others produce 
misleading dynamics when integrate Ordinary Differential Equations (ODEs). 
Typical problems are, for instance, the convergence to spurious steady states, 
changes in the basis of attraction, appearance of spurious bifurcations, etc. The 
essence of such difficulties is that the dynamic of the numerical schemes (seems 
as discrete dynamical systems) is far richer than that of its continuous counter- 
parts. Contrary to the popular belief, drawbacks of this type may no be solved by 
reducing the stepsize of the numerical method. Therefore, it is highly desirable 
the development of numerical integrators that preserve, as much as possible, the 
dynamical properties of the underlaying dynamical system for all step sizes or 



Preprint submitted to a journal 



September 10, 2012 



relative big ones. In this direction, some modest advances has been achieved by a 
number of relative recent integrators of the class of Exponential Methods, which 
are characterized by the explicit use of exponentials to obtain an approximate so- 
lution. An example of such integrators are the High Order Local Linearization 
(HOLL) methods based on Runge-Kutta schemes [4, 5, 14]. 

HOLL integrators are obtained by splitting, at each time step, the solution of 
the original ODE in two parts: the solution of a linear ordinary differential equa- 
tion plus the solution of an auxiliary ODE. The first one is solved by a Local Lin- 
earization (LL) scheme [12, 13] in such a way that A-stability is ensured, while the 
second one can be approximated by any conventional scheme, preferably a high 
order explicit scheme. Originally, HOLL methods were introduced as a flexible 
approach for increasing the order of convergence of the order-2 LL method but, 
in addition, they can be thought as a strategy for constructing high order A-stable 
explicit schemes based on conventional explicit integrators. For this reason, if we 
focus on the conventional integrator involved in a particular HOLL method, then 
it is natural to say that the first one has been locally linearized. In this way, if a 
Runge Kutta scheme is used to approximate the above mentioned auxiliary ODE, 
the resulting HOLL scheme are indistinctly called Local Linearization - Runge 
Kutta (LLRK) scheme or Locally Linearized Runge Kutta (LLRK) scheme. 

In [5], general results on the convergence, stability and dynamical properties 
of the Locally Linearized Runge Kutta schemes were demonstrated. Simulation 
studies carried out in [4, 5, 20] have shown that, for a variety of test equations, 
LLRK schemes of order 3 and 4 preserve much better the stability and dynam- 
ical properties of the actual solutions than their corresponding conventional RK 
schemes. 

However, the accuracy and computational efficiency of the Local Lineariza- 
tion methods have been much less considered up to now, being the dynamical 
properties of such schemes the focus of previous studies and the main reason for 
the development of these methods. The few available results are the following. 
On an identical time partition [5], the LLRK scheme based on the classical order- 
4 RK scheme displays better accuracy than the order-5 RK formula of Dormand 
& Prince [6] in the integration of a variety of ODEs. On different time partitions 
[20], similar results are obtained by an adaptive implementation of the mentioned 
LLRK scheme in comparison with the Matlab code ode45, which provides an 
adaptive implementation of the embedded RK formulas of Dormand & Prince. 
However, this is achieved at expense of additional evaluations of the vector field, 
and with larger overall computational time. With this respect, the main draw- 
back of that adaptive LLRK scheme is the absence of a computationaly efficient 
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strategy based on embedded formulas. 

The purpose of this work is introducing adaptive LLRK schemes based on the 
embedded RK formulas of Dormand & Prince, and evaluating their accuracy and 
computational efficiency in order to study the effect of the local linearization on 
these RK formulas. 

The paper is organized as follows. In the first section, a basic introduction on 

the Local Linearization - Runge Kutta (LLRK) schemes is presented. In section 
2, the embedded Locally Linearized Runge-Kutta formulas are defined, and two 
adaptive implementations of them are described. In the last two sections, the 
results of a variety numerical simulations are presented and discussed respectively. 



2. Notations and Preliminaries 

Let ^ c be an open set. Consider the J-dimensional differential equation 

^ = f(?,x(0), te[to,T] (1) 
x(?o) = xo, (2) 

where xq G ^ is a given initial point, and f : [to, T]x^ — ^ is a differentiable 
function. Lipschitz and smoothness conditions on the function f are assumed in 
order to ensure a unique solution of this equation in ^. 

Let (t) = {tn : n = 0,1, ... ,N} he a time discretization with maximum step- 
size h defined as a sequence of times that satisfy the conditions ?o < < • • • < 
ttf = T and sup(/j„) <h<\, where hn = ?n+i —tniorn = Q,...,N—\. 

n 

For a given (?„,y„), let v„+i = v„ + Ai v„;/j„) be an order-71 approximation 
to solution of the linear ODE 

= B„zi(0 + b„(0, re(?„,?„+i], (3) 
= y« (4) 

at tn+\, and let w^+i = w„ + A2" {tniyfn\hn) be an order- )^ Runge-Kutta scheme 
approximating the solution of the nonlinear ODE 

= q(?n,Zl(?«);^Z2(0)> te{tn,tn+\], (5) 
Z2(?n) = (6) 
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at tn+i, where B„ = fx y„) is ad x d constant matrix, and 

hn{t) = it {tn,yn) {t - tn) + f y«) " B„y„ 

and 

q{tn, y„; 5, <^ ) = f{s, y„ + Ai y„; 5 - + <^ ) - fx(?«, y«)Ai y^; * - 
-f?(?«,y«) -f(?«,y«) 

are J-dimensional vectors. Here fx and ff denote partial derivatives respect to 
X and t, respectivelly. 

Definition 1. ([4, 14, 5]) A Local Linearization - Runge Kutta (LLRK) scheme for 
the ODE (l)-(2) is defined by the recursive expression 

y«+i =y„+Ai(?„,y«;/?„)+A2"(?n,0;/i„) (7) 

for all tn G {t)fj, starting with yo = xq. 

The convergence rate of such numerical integrators is stated in the following 
theorem. 

Theorem 2. ([5]) Let x be the solution of the ODE (l)-(2) with vector field f G 
^^+i([fo,r] X ^,W'), where Y = mm{Yi,Y2}. Then, the LLRK scheme (7) satis- 
fies that 

l|x(wi) -y«+ill = o{h'^) 

for all tn+i G {t)j^ and h small enough. 

According to the Definition 1 , a variety of LLRK schemes can be derived. In 
previous works [4, 5], the Local Linearization scheme based on Pade approxima- 
tions [12, 13] has been used to integrate the linear ODE (3)-(4), whereas the so 
called four order classical Runge-Kutta scheme [2] has been applied to integrate 
the nonlinear ODE (5)-(6). This yields the order-4 LLRK scheme 

hn 

y«+i =y« + U4 + ^(2k2 + 2k3+k4), (8) 

where 

u, = L(P6,6(2-'^^D„c,/z„))2'^r 



4 



and 



k^- = f + cy/z„ , y„ + + c;/z„ky_ 1 ) - f (?„ , y„ ) - fx y„ ) - ff (r„ , y„ ) c;/z„ , 

with ki = and c = [ ^ ^ 1 ] . Here, Pp,^(-) denotes the (j!?,^)-Pade ap- 
proximation for exponential matrices [16], and Kj the smallest integer number 
such that ||2^'^>D„Cj/j„|| < The matrices D„, L and r are defined as 



fx(?n,y«) f?(?n,yn) f(?n,y«) 

1 





+2)x{d+2) 



L, = [id 0^x2 ] and = [ Oix(rf+i) 1 ] for of non-autonomous ODEs; and 
as 



fx(y«) f(yn) 




,(rf+l)x(d+l) 



L = [ 0^x1 ] and rT = [ Oi^d 1 ] for autonomous equations. 

On an identical time partition [5], LLRK formula (8) displays better accuracy 
than the order- 5 RK formula of Dormand & Prince in the integration of a variety 
of ODEs. On different time partitions [20], similar results are obtained by an 
adaptive implementation of LLRK formula (8) in comparison with the Matlab 
ode45 code, which provides an adaptive implementation of the embedded RK 
formulas of Dormand & Prince. However, this is achieved at expense of additional 
evaluations of the vector field f, and with larger overall computational time. This 
cost can be sensitively reduced by using the (2, 2)-Pade approximations instead of 
the order (6,6) one used in formula (8), preserving the order of convergence and 
without significant lost of accuracy [20]. 



3. Numerical scheme 

3. 1 . Embedded Locally Linearized Runge-Kutta formulas 

In view of the Definition 1 , new integration formulas can be obtained as fol- 
lows. Similarly to the LLRK scheme (8), the Local Linearization scheme based 
on Pade approximations [12, 13] is used for integrating the linear ODE (3)-(4) but, 
instead of the classical order-4 RK scheme, the embedded Runge-Kutta formulas 
of Dormand & Prince [6] is now applied to integrate the nonlinear ODE (5)-(6). 
This yields the embedded Locally Linearized Runge-Kutta formulas 
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s s 

yn+1 = yn + u,- + hnY,b jk j and y„+i = y„ + + /z„ bjk j, (9) 
where 5 = 6 is the number of the stages, 

and 

i=l 

with ki = and Runge-Kutta coefficients ay bj, bj and cj defined in the Table 
1. Here, Pp,g(-) denotes the (p,^)-Pade approximation for exponential matrices 
with p + q> 4. The number Kj and the matrices D„, L and r are defined as in the 
previous section. 

According with Lemma 4. 1 in [15], if the (/?, q)-Fad6 approximation for expo- 
nential matrices is used, then the order 71 in Theorem lis yi = p + q. Thus, since 
p + q > 4, the order of convergence of LLRK formulas yn+i and y„_|_i are 5 and 
4, respectively. 



5 5 

3 9 

10 40 40 

4 44 _56 32 

5 45 "Ts T 
19372 25360 64448 212 

9 6561 ^ 2187 6561 ^729 

y017 355 46732 49 5103 



3168 33 5247 176 18656 

35 500 125 2187 n_ 

384 1113 192 ~6784 84 

35 500 125 2187 n_ 

384 1113 192 ~6784 84 

5179 7571 393 92097 187 _l_ 

57600 16695 640 ~ 339200 2100 40 



Figure 1: Coefficients tableau for the embedded formulas. 
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3.2. Adaptive strategy 

In order to write a code that automatically adjust the stepsizes for achiving a 
prescribed tolerance of the local error at each step, an adequate adaptive strategy 
is necessary. At glance, the automatic stepsize control for the embedded RK for- 
mulas of Dormand & Prince seems to fit well for the embedded LLRK formulas 
(9). In what follows, the adaptive strategy of the Matlab code ode45 for these 
formulas is described. 

Once the values for the relative and absolute tolerances RTol and ATol, and 
for the maximum and minimum stepsizes h^^x and /imin are set, the basic steps of 
the algorithm are: 

1. Estimation of the initial stepsize 



min{/i, 



max{/?min,A}} 



where 




with 




and tr = ^f^. Initialize fail = 0. 

2. Evaluation of the embedded formula (9) 

3. Estimation of the error 




i=l,...,d 



4. Estimation of a new stepsize 



■new — 



min{/z, 



•max 5 



max{/7min,A}} 



where 




0.8-(|^)^/^-/z if error < RTol 

max|0.1,0.8- (^)^^^l if error > RTol and fail = 

I- ' ^ error 'J •' 

0.5 ■ h if error > RTol and fail — 1 
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5. Validation of yn+i- if error < RTol, then accept as an approximation 
to X at tn+i = tn + h. Otherwise, return to 2 with hn = knew and fail = 1 . 

6. Control of the final step: if tn + h = T , stop. If tn + h + knew > T, then 
redefine hng^ = T — {tn + h) . 

1. Return to 2 with n = n-\-\,hn = h„ew, and fail = 0. 

3.3. Continuous formula 

Continuous formulas of RK methods are usually defined for computing the 
solutions on a dense set of time instants with minimum computational cost. Typi- 
cally [7], they are constructing by means of a polynomial interpolation of the RK 
formulas between two consecutive times ?„,?„+i G (f)^. 

By a simple combination of the LLRK formulas (9) with the continuous for- 
mula of the Dormand & Prince RK method [7] for (5)-(6), a continuous 5-stage 
LLRK formula can be defined as 

s 

y{t„ + eK) = yn+u{eh„) +h„ £ bj{t„ + eh„)kj , o < e < i, (lO) 

for all tn,tn+i e (O/i' where 
is a J-dimensional vector, and 

;=1 

is a polynomial with coefficients a, j. Here, the function k^, the matrices D„, L 
and r, and the number Kj are defined as in (9), as well as the (/7,^)-Pade approxi- 
mation Pp g. The coefficients atj, defined in Table 1, coincides with those of the 
continuous RK formula implemented in the Matlab code ode45. 

3.4. LLRK45code 

This subsection describes a Matlab2007b(32bits) implementation of the adap- 
tive scheme described above, which will be denoted as LLRK45 code. 

In order to make a fair comparison with the Matlab code ode45, the LLRK45 
code is an exact copy of the ode45 one with the exception of the program lines 
corresponding to the embedded and continuous formulas of Dormand and Prince, 
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j/i 


1 


2 


3 


4 


1 


1 


-183/64 


37/12 


-145/128 


2 














3 





1500/371 


-1000/159 


1000/371 


4 





-125/32 


125/12 


-375/64 


5 





9477/3392 


-729/106 


25515/6784 


6 





-11/7 


11/3 


-55/28 


7 





3/2 


-4 


5/2 



Table 1: Values of the coefficient a, j involved in the continuous LLRK formula (10) . 

which are replaced by the formulas (9) and (10) respectively. We recall that the 
code ode45 implements the adaptive strategy of the subsection 3.2 for the embed- 
ded RK formulas of Dormand & Prince, which is considered for many authors the 
most recommendable code to be applied as first try for most problems [17]. 

Note that, the embedded LLRK formulas (9) require the computation of six 
Fade approximations P^.^ at each integration step, which increases the computa- 
tional cost of the original embedded RK formulas. Nevertheless, this number of 
Fade approximations can be reduced by taking in to account that: a) (Pp.^(2^'^D„cj/j„))^ 
gives an approximation to exponential matrix ^nCjh„.^ j-j^g property of 

the exponential operator. Indeed, this can be carried out in two steps: 

1. approximating e°"''"/^° by the matrix Mi/90 = (Pp,^(2"'^D„/7,,/90))2'', where 
K is the smallest integer number such that ||2^'^D„/z„/90|| < |; and 

2. the successive computation of the matrices 



M2/90 


= M1/90M1/90 


M4/90 = 


M2/90M2/90 


M8/90 


= M4/90M4/90 


M16/90 = 


M8/90M8/90 


M32/9O 


= M16/90M16/90 


M80/90 = 


= M32/90M16/90M32/90 


Ml/10 


= M8/90M1/90 


Mi/5 = 


Mi/ioMi/10 


M2/5 


= M1/5M1/5 


M4/5 = 


M2/5M2/5 


M3/10 


= Mi/ioMi/5 


Ml = 


M4/5M1/5. 



Consequently, the matrix Mf^. corresponding to each RK coefficient Cj pro- 
vides an approximation to e^"*^'^", for all j = 1, ..,6. In this way, at each integra- 
tion step, the code LLRK45 performs six evaluation of f (same than the ode45 
code), one Jacobian matrix and one matrix exponential. 
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The matrix Mi/90 is computed by means of the function "expmf", which pro- 
vides a C++ implementation of the classical (p, ^)-Pade approximations algorithm 
for exponential matrices with scaling and squaring strategy [16], and p = q = 3. 



4. Numerical simulations 

In this section, the performance of the LLRK45 and ode45 codes is compared 
by means of numerical simulations. To do so, a variety of ODEs and simulation 
types were selected. For all of them, the Relative Error 



RE = max 

i=l,...,d;tj<E{t)i, 



x'itj) 



(11) 



between the "exact" solution x and its approximation y is evaluated. 

We recall that the code ode45 is less efficient at crude tolerance (e.g., at the 
default tolerance RTol = 10^ and ATol = 10^^) [17], therefore we fix the mild 
tolerance RTol = 10^ and ATol = 10^^ for this code in all simulations. The 
Matlab code odelS* with refined tolerance RTol = 10^^ md ATol = 10~^^ is 
used to compute the "exact" solution x in all simulations. 

4.1. Test Examples 

The first four examples have the semi-lineal form 

dx 

5^=Ax + g(x), (12) 

where A is a square matrix and g is a function of x. The vector field of the first two 
examples have Jacobians with eigenvalues on or near the imaginary axis, which 
made these equations difficult to be integrated by conventional schemes [17]. The 
other two are also hard for conventional explicit schemes since they are examples 
of stiff equations [17]. Example 6 has an additional complexity for a number of 
integrators that do not update the Jacobians of the vector field at each integration 
step [17, 11]: the Jacobian of the linear term has positive eigenvalues, which 
results a problem for the integration in a neighborhood of the stable equilibrium 
point X = 1. 

Example 3. Periodic linear [5 ] 

|^A(x + 2), 
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with 

—I 

x{to) = (-2.5,-1.5) and [to,T] = [0,4n]. 
Example 4. Periodic linear plus nonlinear term [5 ] 

^=A(x + 2)+0.1x2, 
dt 

where the matrix A is defined as in the previous example, x(?o) = (1,1). cind 

[tQj] = [OATi]. 

Examples. Stiff linear [5] 

dx 

— = -100H(x + l), 
dt 

where H is the 12 -dimensional Hilbert matrix (with conditioned number 1.69 x 
IQi^j, x'(?o) = 1, / = 1 . . . 12, and [?o, T] = [0, 1]. 

Example 6. Stiff linear plus nonlinear term [5 ] 
dx 

— = 100H(x-l) + 100(x-l)2-60(x3-l), 
dt 

where H is the 12-dimensional Hilbert matrix, x'(?o) = —0.5, i = 1 ... 12, and 

[?o,r] = [o,i]. 

The following examples are well known nonlinear test equations. This include 
highly oscillatory, non stiff and mild stiff equations. 

Example 7. Fermi-Pasta-Ulam equation defined by the Hamiltonian system [9] 

J 3 2 3 3 

^(P, q) = ;^ E(pLi +p|-) + X I^(^2, - q2/-i)' + i;(q2,+i - q2,)^ 

with w = 50, initial conditions 1, 1, 1/w, 1 for the four first variables and zero for 
the remainder eight, and [to^T] = [0, 15]. 
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Example 8. Brusselator equation [7]: 

— = 3xi-XiX2, 

at 

where (xi(to),X2(to)) = (1-5,3) and [to^T] = [0,20]. 

Example 9. Rigid body equation [7]: 

dxi 
dt 

dX2 

dt 

= —0.51x1X2 

dt 

with {xi{to),X2{to),X3{to)) = (0, 1, 1) over [to,T] = [0, 12]. 

Example 10. Chemical reaction [17]: 

dx\ 
dt 

dX2 

dt 
dx3 

dt 

dX4 



,^ X2X3 
dt 



,^ -X1X3 
dt 



1.3fe-JCi) + 10400A:(xi)jC2 

dt 

-- 1880(jC4-JC2(l+^(xi))) 

: 1752-269x3 + 267x1 

dt 



-- 0.1 + 320x2-321x4 
dt 



(20 7 '^"" ) 

where k{xi) = e ' . With initial condition (50,0,600,0.1) over [to,T] ~ 
[0, 1], this is mild stijf equation. 

Example 11. Van der Pol equation [8]: 

dx\ 

dX2 2\ 

— = e(l-X2)Xi+X2 
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Example 


Code 


RTol 


ATol 


Time 
steps 


Relative 
Error 


PerLin 


ode45 


10-* 


10-^^ 


598 


2.6 X 10-^ 


LLRK45 






598 


2.0 X 10-" 


PeiNoLin 


ode45 


lO-*" 


io-« 


411 


2.6 X 10-^ 


LLRK45 






411 


6.9 X lO-'' 


StiffLin 


ode45 


10-" 


io-« 


78 


1.1 X 10-" 


LLRK45 






78 


2.7 X IQ-'^ 


StiffNoLin 


ode45 


lO-*" 


10-^^ 


133 


1.5 X 10-5 


LLRK45 






133 


6.8 X 10-" 



Table 2: Relative error of the order-5 formula of each code when integrate the semilinear examples 
over identical time partition. 

with {xi{tQ),X2{to)) = (2,0). With e = 1 and £ = 10^, this is a non stiff anda mild 
stiff equation on the intervals [to,T] = [0,20] and [tQ,T] = [0,300], respectively. 

As illustration, Figure 2 shows the first component of the solution of each ex- 
ample, which will be consecutively named as PerLin, PerNoLin, StiffLin, StiffNo- 
Lin, fpu, bruss, rigid, chm, vdpl and vdplOO. 

4.2. Simulation A: Integration over same time partition 

This simulation is designed to compare the accuracy of the order-5 formulas 
of the codes LLRK45 and ode45 over identical time partition. First, the ode45 
code integrates all the examples with the tolerances RTol = 10^ and ATol ~ 
10^^. This defined, for each example, a time partition {t)h over which the order-5 
formula of the LLRK45 code is evaluated as well. That is, the formula 

y«+i = Yn + + hn b jk j 
7=1 

with = L(P3,3(2"'^/D„c/7„))2''^r for the LLRK45 code. Tables 2 and 3 present, 
respectively, the Relative Error (11) of each code in the integration of the four 
semilinear and six nonlinear examples defined above. The number of accepted 
time steps is shown as well. 

4.3. Simulation B: Integration with same tolerance 

This simulation is designed to compare the performance the codes LLRK45 
and ode45 with the same mild tolerances RTol = 10^^ and ATol = 10^^. As 
a difference with Simulation A, here each codes use a different time partition 
defined by the adaptive strategy. 
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Figure 2: Path of the first components of the solution in each example. 
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Example 


Code 


RTol 


ATol 


Time 


Relative 


fpu 


ode45 


10-^ 


10-" 


4474 


8.1 


LLRK45 






4474 


2.9 X 10-^' 


rigid 


ode45 


10-'' 


10-'-- 


66 


7.5 X 10-^ 


LLRK45 






66 


4.0 X 10-" 


chm 


ode45 


10-^ 


10-" 


723 


1.1 X 10-" 


LLRK45 






723 


2.5 X 10-' 


brass 


ode45 


lO-*" 


10-" 


148 


8.7 X 10-" 


LLRK45 






148 


3.5 X 10-' 


vdpl 


ode45 




10-" 


204 


2.8 X lO-"* 


LLRK45 






204 


1.5 X 10-' 


vdplOO 


ode45 


10-^ 


10-" 


17516 


0.42 


LLRK45 






17516 


9.2 X 10-^ 



Table 3: Relative error of the order-5 formula of each code when integrate the nonlinear examples 
over identical time partition. 



Example 


Code 


Time 
steps 


Failed 
steps 


/ 


/v 


Expm 


Time 


Relative 
EiTor 


PerLin 


ode45 


598 





3589 








1 


2.6 X 10--^ 




LLRK45 


14 





85 


14 


14 


0.08 


3.0 X 10-" 


PerNoLin 


ode45 


411 





2467 








1 


2.6 X K)-' 




LLRK45 


137 





823 


137 


137 


0.62 


3.6 X 10-" 


StiffLin 


ode45 


78 


1 


475 








1 


1.1 X 10-" 




LLRK45 


14 





85 


14 


14 


0.34 


2.3 X 10-'^ 


StiftNoLin 


ode45 


133 


5 


829 








1 


1.5 X 10-' 




LLRK45 


43 





259 


43 


43 


0.53 


1.6 X 10-" 



Table 4: Code performance in the integration of the semilinear examples with the same tolerances. 

Tables 4 and 5 summarize the results of each code in the integration of each 
example. The column "Time" in these tables presents the relative overall time 
of each numerical scheme with respect to that of the ode45 code on the whole 
interval [^o,?"]- This overall time ratio works as simple indicator to compare the 
total computational cost of each code. In addition, the tables show the number of 
accepted and failed steps, the number of evaluations of f and f^, and the number 
of exponential matrices computed in the integration of each example. 

4.4. Simulation C: Integration with similar accuracy 

In this type of simulation, the tolerances RTol and ATol of the LLRK45 code is 
changed until its relative error in the integration of each example achieves similar 
value to that corresponding to the code ode45. Tables 6 and 7 summarize the 
performance of each code in the integration of each example. 
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Example 


Code 


Time 
steps 


Failed 
steps 


/ 


/.■ 


Expm 


Time 


Relative 
Error 


fpu 


ode45 


4474 


60 


27205 








1 


8.1 


LLRK45 


1496 


125 


9727 


1621 


1621 


0.67 


2.0 X 10"^ 


rigid 


ode45 


66 


4 


421 








1 


7.5 X 10-5 


LLRK45 


53 


5 


349 


58 


58 


1.55 


8.6 X 10-" 


chm 


ode45 


723 


16 


4435 








1 


1.2 X IQ-" 


LLRK45 


357 


2 


2155 


359 


359 


1.05 


9.2 X IQ-' 


brass 


ode45 


148 


13 


967 








1 


8.7 X lO-*" 


LLRK45 


105 


14 


715 


119 


119 


1.47 


5.4 X 10-" 


vdpl 


ode45 


204 


32 


1417 








1 


2.8 X lO-"* 


LLRK45 


162 


38 


1201 


200 


200 


1.45 


5.8 X 10-5 


vdplOO 


ode45 


17516 


1540 


114337 








1 


0.42 


LLRK45 


7893 


19 


47473 


7912 


7912 


0.69 


2.1 X 10"^ 



Table 5: Code performance in the integration of the nonUnear examples with the same tolerances. 



Example 


Code 


RTol 


ATol 


Time 

steps 


Failed 

steps 


/ 




Expm 


Time 


Relative 
Error 


PerLin 


ode45 


10"" 


10-*^ 


598 





3589 








1 


2.6 X 10-^ 




LLRK45 


0.1 


lO-'' 


13 





79 


13 


13 


0.06 


2.0 X K)-" 


PerNoLin 


ode45 


10-" 


10-y 


411 





2467 








1 


2.6 X 10-5 




LLRK45 


7 X 10"" 


7 X 10-^ 


95 





571 


95 


95 


0.44 


1.5 X 10-5 


StifELin 


ode45 


10-" 


10-" 


78 


1 


475 








1 


1.1 X K)-" 




LLRK45 


0.1 


10-* 


13 





79 


13 


13 


0.31 


2.6 X 10-" 


StiffNoLin 


ode45 


10-" 


10-" 


133 


5 


829 








1 


1.5 X 10-5 




LLRK45 


4 X 10-5 


4x 10-" 


28 





169 


28 


28 


0.35 


1.1 X 10-5 



Table 6: Code performance in the integration of the semilinear examples with similar accuracy. 



Example 


Code 


RTol 


ATol 


Time 
steps 


Failed 
steps 


/ 


/>- 


Expm 


Time 


Relative 
Error 


fpu 


ode45 


10-" 


10-" 


4474 


60 


27205 








1 


8.1 


LLRK45 


10-'' 


10-' 


567 


39 


3637 


606 


606 


0.26 


5.2 


rigid 


ode45 


10-" 


10-" 


66 


4 


421 








1 


7.5 X 10-5 


LLRK45 


2.2 X 10-5 


2.2 X lO-*' 


30 


2 


193 


32 


32 


0.93 


3.3 X 10-5 


chm 


ode45 


10"" 


10-" 


723 


16 


4435 








1 


1.1 X K)-" 


LLRK45 


1.4 X 10"" 


1.4 X 10-^ 


341 


2 


2059 


343 


343 


1 


1.0 X 10-" 


brass 


ode45 


10-" 


10-" 


148 


13 


967 








1 


8.7 X 10-" 


LLRK45 


1.2 X 10-" 


1.2 X 10-" 


101 


13 


685 


114 


114 


1.27 


8.1 X 10-" 


vdpl 


ode45 


10-" 


10-" 


204 


32 


1417 








1 


2.8 X 10-"* 


LLRK45 


3.2 X K)-" 


3.2 X 10-" 


128 


26 


925 


154 


154 


1.12 


2.4 X lO-"* 


vdplOO 


ode45 


10-" 


10-" 


17516 


1540 


114337 








1 


0.42 


LLRK45 


3.9 X 10-5 


3.9 X 10-** 


5026 


731 


34543 


5757 


5757 


0.50 


0.26 



Table 7: Code performance in the integration of the nonlinear examples with similar accuracy. 
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Example 


Code 


Time 
steps 


Dense 
OutPut 


Relative 
Error 


PerLin 


ode45 


598 


2393 


2.6 X 10"- 


LLRK45 


14 


57 


3.0 X IQ-'' 


PerNoLin 


ode45 


411 


1169 


3.0 X ur^ 


LLRK45 


137 


293 


8.7 X ur' 


StiffLin 


ode45 


78 


313 


1.1 X 10"^ 


LLRK45 


14 


57 


3.4 X 10-'^ 


StiffNoLin 


ode45 


133 


533 


3.0 X IQ-^ 


LLRK45 


43 


173 


2.9 X 10"' 



Table 8: Relative error of the continuous formulas of the codes over their dense output after 
integrating the semilinear examples. 



Example 


Code 


Time 
steps 


Dense 
OutPut 


Relative 
Error 


fpu 


ode45 


4474 


17897 


19.0 


LLRK45 


1496 


5985 


2.8 X IQ-^ 


rigid 


ode45 


66 


265 


3.4 X Kr"* 


LLRK45 


53 


213 


1.7 X ur" 


chm 


ode45 


723 


2893 


1.1 X lO-*" 


LLRK45 


357 


1429 


9.2 X ur' 


bruss 


ode45 


148 


593 


1.0 X 10-=' 


LLRK45 


105 


421 


2.4 X 10-' 


vdpl 


ode45 


204 


817 


6.9 X lO-'' 


LLRK45 


162 


649 


2.3 X ur" 


vdpl 00 


ode45 


17516 


70065 


0.47 


LLRK45 


7893 


31573 


4.1 X 10-^ 



Table 9: Relative error of the continuous formulas of the codes over their dense output after 
integrating the nonlinear examples. 

4.5. Simulation D: Evaluation of the dense output 

This simulation is designed to compare the acuracy of the continuous formulas 
of the codes LLRK45 and ode45 over their dense output. For this, both codes are 
applied to each example with the same tolerances RTol = 10^ and ATol = 10^ 
but, the relative error of each code is now computed on its respective dense output 
instead on the time partition {t)h defined by the adaptive strategy. Tables 8 and 9 
present these relative errors. The number of accepted time steps and dense output 
times are also shown. 

5. Discussion 

The results of the previous section show the following. On the same time 
partition, the embedded LLRK formulas are significantly much accurate than the 
classical embedded RK formulas of Dormand & Prince. With identical tolerances 
and adaptive strategy, LLRK45 code is much accurate than the ode45 code and 
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requires much less time steps for integrating the whole intervals. The overall time 
of the adaptive LLRK45 code is, in general, lower or similar than the overall time 
of the ode45 code. The accuracy of the dense output of the LLRK code is, in 
general, higher than the accuracy of the ode45 code. 

We recall that, in order to study the effect of the local linearization on the 
conventional RK scheme of Dormand and Prince, the LLRK code considered in 
this work is an exact copy of the code ode45 with the exception of the program 
lines corresponding to the embedded and continuous formulas. In this way, the 
code does not include a number of necessary modifications that might improve 
the performance of the local linearized integrator. Some of they are the following: 

• the initial h atto, which can be estimated by means the exact second deriva- 
tive of the solution x with no extra cost [20]; 

• online smoothness and stiffness control for estimating the new h at each step 
(as, e.g., in [8]); 

• the automatic detection of constant Jacobian matrix (as in [17]); 

• option for using exact, numerical or automatic Jacobian matrices (as in [17, 
18, 1]); 

• faster algoritms to compute the Pade approximation to exponential matrix 
(as, e.g., in [10]) 

• a parallel implementation of matrix multiplications involved in the expo- 
nential matrix evaluations for taking advantage of the multi core technology 
available in the current microprocessors; 

• increase the number of times of the dense outputs: a) up to twelve per each 
pair of consecutive times of the partition {t)h with no extra computation of 
exponential matrices; or b) up to ninety with some few extra matrix multi- 
plications; 

• a new continuous formula that replace the current one based on the continu- 
ous RK formula by other based on a polynomial interpolation of the LLRK 
formula itself on the whole interval (i.e, derived from the standard way of 
constructing continuous RK formulas as in [7]); and 

• change of hmax^ which seems to be too short for semilinear equations. 
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6. Conclusions 



The local linearization of the embedded Runge-Kutta formulas of Dormand 
& Prince produces a significative improvement of the accuracy of these formulas, 
which implies a substantial reduction of the number of time steps and, conse- 
quently, a sensitive reduction of the overall computation cost of their adaptive 
implementation. 
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